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ABSTRACT 

The  lecture  addresses  topics  in  multivariate  approximation  which  have  caught  the 
author’s  interest  in  the  last  ten  years.  These  include:  the  approximation  by  functions 
with  fewer  variables,  correct  points  for  polynomial  interpolation,  the  B(ernstein,-ezier,- 
arycentric)-form  for  polynomials  and  its  use  in  understanding  smooth  pp  functions,  approx¬ 
imation  order  from  spaces  of  pp  functions,  multivariate  B-splines,  and  surface  generation 
by  subdivision. 
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SIGNIFICANCE  AND  EXPLANATION 


This  lecture  was  given  at  the  IMA  'SIAM  conference  on  “State  of  the  Art  in  Numerical 
Analysis’’  held  April  14-18.  1986,  in  Birmingham,  England.  The  lecture  reviews  develop¬ 
ments  in  Multivariate  Approximation  in  the  last  ten  years.  The  selection  of  topics  is  quite 
subjective;  it  reflects  entirely  the  author's  research  experience  during  that  time. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary  lies 
with  MRC,  and  not  with  the  author  of  this  report. 


MULTIVARIATE  APPROXIMATION 
Carl  de  Boor 


1.  The  set-up 

The  talk  concerns  the  approximation  of 

/  :  G  C  Ud  — *  IR, 

i.e.,  of  some  real-valued  function  /  defined  on  some  domain  G  in  d-dimensional  space. 

While  my  first  publication  dealt  with  such  multivariate  (actually,  bivariate)  approxi¬ 
mation,  I  have  concerned  myself  seriously  with  multivariate  approximation  only  in  the  last 
ten  years.  This  talk  reflects  some  of  the  experiences  I  have  had  during  that  time. 

The  approximating  functions  are  typically  polynomials  or  piecewise  polynomials,  and 
just  how  one  describes  them  will  have  an  effect  on  one’s  work  with  them.  Papers  on 
multivariate  approximation  often  sink  under  the  burden  of  cumbersome  notation.  As  a 
preventive  measure  against  such  a  sad  fate,  I  shall  follow  the  “default”  convention  whereby 
symbols  are  left  out  if  they  can  reasonably  be  guessed  from  the  context.  For  example,  if 
x  =  (x(l), . . . ,  x(d))  has  just  been  declared  to  be  a  point  in  IRrf,  I  will  feel  free  to  write 

d 

^x(i)  instead  of  ]jT^x(t). 

i  t=l 

For  polynomials,  multi-index  notation  is  standard.  With  x  =  (x(l), . . .  ,x(d))  the 
generic  point  in  IR**,  one  uses  the  abbreviation 

x°  :=  I ]  x(Oo(0,  x  €  IRd  ,  a  €  TLd+. 

i 

The  function  x  >-*  xa  is  a  monomial  of  degree  a,  or,  of  (total)  degree 

|q|  :=  5^a(i) 

1 

if  only  the  exponent  sum  matters.  More  generally,  a  polynomial  of  degree  <  a  is,  by 
definition,  any  function  of  the  form 

x  ~  ^  xpc(3), 

&<.Q 
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with  real  coefficients  c(0).  The  collection  of  all  such  polynomials  is  denoted  by 


7Ta  =  7T0(IRd). 

the  collection  of  all  polynomials  of  total  degree  at  most  k  by 

7Tfc  =  7T*(IRd), 

and  the  collection  of  all  polynomials,  of  whatever  degree,  by 

7T  =  7r(IRd). 

Many  expressions  simplify  if  one  makes  use  of  the  normalized  power  function,  i.e., 
the  function 

0 a:x~x°/a\  :=  J]*(t)-(,)/a(i)!. 

t 

For  example,  with  a,  £,  v, . . . ,  f  6  ZZ^ ,  the  Multinomial  Theorem  takes  the  simple  form 

I* y  (ii) 

£  1 /  .  -f  £  —  0* 

At  times,  it  pays  to  give  up  on  the  power  form  altogether.  In  some  contexts,  it  is  very 
convenient  to  describe  polynomials  in  terms  of  the  particular  homogenous  polynomials 

(K,  •)  :  x  —  <y,x), 

with 

<Vi*>  5^»(*)x(0 

X 

and  Y  a  finite  subset  of  IR^.  I  will  make  use  of  this  form  in  Section  3.  A  particular  instance 
is  the  B(ernstein-ezier)-form,  which  is  the  form  of  choice  when  dealing  with  piecewise 
polynomials  on  a  triangulation  (see  Section  4). 

2.  Approximation  by  functions  of  fewer  variables 

The  simplest  approach  to  multivariate  approximation  uses  tensor  products,  i.e., 
linear  combinations  of  functions  of  the  form 

t — n 
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each  gt  being  a  univariate  function.  This  neatly  avoids  dealing  with  the  realities  of  multi¬ 
variate  functions,  but  its  effectiveness  depends  on  having  the  information  about  /  corre¬ 
spondingly  available  in  (cartesian)  product  form,  e.g.,  as  function  values  on  a  rectangular 
grid  parallel  to  the  coordinate  axes.  The  recent  book  by  Light  and  Cheney  (1986)  provides 
up-to-date  material  on  this  practically  very  important  choice  of  approximating  function 
and  certain  ready  extensions. 

The  book  also  deals  with  the  situation  when  the  information  about  /  is  not  in  product 
form.  In  that  case,  tensor  product  approximants  still  are  attractive  since  they  are  composed 


of  univariate  functions.  The  general  question  of  how  to  approximate  a  multivariate  function 
by  functions  of  fewer  variables  has  received  much  attention.  A  reference  with  a  Numerical 
Analysis  slant  is  Golomb  (1959). 

The  most  remarkable  result  along  this  line  is 

Kolmogorov’s  Theorem  (Kolmogorov,  1957) 

3A  6)0,  l]d,  3$  C  Lipa|0, 1]  n  strictly  monotone,  =  2d  +  1,  such  that 
V/€C[ 0,l]d  3geC\0,d] 

/(*)  =  ^  *('52  *(*)*>(*(*)))• 

i 


(Here  and  below, 


I 


#A  :=  the  number  of  elements  in  A.) 


The  theorem  claims  the  existence  of  a  set  if  of  2 d+ 1  ‘universal’  maps  tji  :  [0,  l]d  — ♦  [0,  d] 
so  that,  for  each  continuous  function  /  on  the  unit  cube  jO,  l]d,  a  continuous  function  g  on 
the  interval  [0,  d)  can  be  found  for  which 

/  =  £<7°0. 

Moreover,  each  function  ip  6  'k  is  of  the  form 

d 

ii’(i)  :=  ^  A(i)v?(j(t)), 


mm 


•./'A  VrV/.'.V^.vyv.V  .■ 'V1.-  V v" .-'vV  '■  V-'-.' >V.\-.\vV-V-\vV- V- 
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with  each  of  the  2d  +  1  functions  <p  6  $  strictly  increasing  and  in  LipQ  0,  lj  for  some 
positive  q,  and  A  some  d-vector  with  positive  entries  all  bounded  by  1. 

The  work  of  Kolmogorov  and  his  pupil  Arnol'd  which  culminated  in  this  theorem 
was  motivated  by  Hilbert’s  Thirteenth  Problem  which  contained  (implicitly)  the  conjec¬ 
ture  that  not  all  continuous  functions  of  three  variables  could  be  written  as  superpositions 
of  continuous  functions  of  two  variables.  The  version  quoted  here  reflects  further  sim¬ 
plifications,  chiefly  by  Lorentz  (1962).  For  a  proof  and  further  discussion,  see  (Lorentz, 
1966;pp.l68ff). 

Practical  use  of  Kolmogorov’s  Theorem  seems  elusive  since  the  ’universal’  functions 
¥3  c  $  have  a  fractal  ’derivative’  (see  Section  8)  and  g  need  not  be  smooth  even  if  /  is 
smooth.  But  it  remains  a  challenge  to  develop  a  practical  Approximation  Theory  which  can 
handle  approximating  functions  of  the  form  (2.1).  In  any  case,  it  suggests  a  nontraditional 
form  of  approximation  which  is  motivated  by  computational  or  algorithmic  simplicity. 
Perhaps  we  have  been  too  accepting  of  traditional  approximation  techniques  in  which  we 
choose  the  approximating  family  according  to  linear  degrees  of  freedom.  Perhaps  we  should 
consider  instead  approximating  families  which  are  classified  by  the  number  of  floating-point 
operations  required  for  their  evaluation.  Approximation  Theory  as  it  now  exists  has  little 
to  offer  in  this  direction,  but  Computer  Science  may  have  something  to  teach  us. 

A  special  case  has  had  much  exposure  in  the  times  before  electronic  computers,  viz., 
the  approximation  by  nomographic  functions.  These  are  functions  of  two  variables  of 
the  specific  form 

IR2  -  -  IR  :  x  ►->  <?(y?(x(l))  -t-  t/>(i(2))). 

For  recent  algorithmic  work,  see  von  Golitschek  (1984). 

3.  Loss  of  Haar 

The  traditional  approaches  to  approximation  all  start  with  polynomials,  and  so  will  1. 
Perhaps  the  greatest  change  when  going  to  the  multivariate  set-up  is  the  loss  of  the  Haar 
property. 
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To  recall,  interpolation  from  a  linear  space  5  of  functions  on  some  G  C  1R  at  a  point 
set  T  C  G  can  be  viewed  as  the  task  of  inverting  the  restriction  map 

5  — *  IR7^  :  p  — *  Pit- 

We  are  to  (re)construct  some  element  p  €  S  from  its  prescribed  values  p  ?  at  the  points  in 
T.  In  these  terms,  S  has  the  Haar  property  if  every  T  C  G  with  4?T  -  dimS  is  correct 
for  5,  i.e.,  is  such  that 

5  — *  IR7  :  p  p\j  is  1-1  and  onto. 

In  other  words,  we  can  interpolate,  and  uniquely  so,  from  5  to  any  function  values  given 
on  T. 

If  d  —  1.  then  ■n^  has  the  Haar  property  (for  any  G  with  more  than  k  points),  but 
Mairhuber’s  switching  yard  argument  (cf.,  e.g.,  the  cover  of  Lorentz  (1966))  shows  that 
this  property  cannot  hold  for  any  S  of  dimension  >  1  on  a  multidimensional  set  G. 

For  the  case  of  polynomials,  it  is  possible  to  identify  various  point  sets  T  C  IRd  which 
are  correct  for  ilk-  A  particularly  nice  example  is  provided  by  Chung  and  Yao  (1977)  who 
prove  the  following:  Suppose  that  V  is  a  finite  subset  of  IRrf\0,  and  that  V  'JO  is  in  general 
position,  which  means  that  7r j  has  the  Haar  property  on  V  j  0.  This  implies  that  for 
every  W  c  V  with  #VF  =  d,  xw  6  IR^  is  defined  uniquely  by  the  equations 

1  +  (w,xw)  =  Or  tt'  G  W, 

since  these  state  that  the  linear  polynomial  1  (-.xw)  is  to  vanish  on  W  and  take  the 

value  1  at  0.  Further,  since  1  +  (-,xw)  already  vanishes  on  the  d  points  in  W,  it  cannot 
vanish  anywhere  on  V'\IF,  by  the  Haar  property.  It  follows  that  the  functions 

TT  1  —  <v.x) 

Lw  ■  i  —  II  v - 

1  -  v.  Xw 

<  6  V  VV 

are  well-defined,  are  made  up  of  -  d  linear  factors,  and  satisfy  iw(iw)  -  (since, 

for  W  *  W.  at  least  one  r  (  \  VV  is  in  W\  making  the  corresponding  linear  factor  zero  at 


xw').  hence  in  particular  xw<  ^  xw  for  Vi"  ^  W',  while  #T  =  (#dV)  =  dim  7r#v-d(IRrf). 
This  proves 

Theorem  (Chung  and  Yao,  1977)  T  :=  {%}  is  correct  for  7r #v-d- 

The  following  result  has  a  different  flavor: 

Theorem  (Hakopian,  1983)  If  T  C  7Ld+  ‘contains  its  shadow’,  i.e.,  0  <  a  G  T 
implies  0  €  T,  then  T  is  correct  for  span^JQ^ 

A  totally  different  approach  to  the  correctness  problem  has  been  taken  by  Kergin 
(1978).  Kergin  is  interested  in  extending  to  a  multivariate  setting  H.  Whitney’s  (1957) 
characterization  of  functions  on  some  subset  T  of  IR  which  have  extensions  to  a  smooth 
function  on  all  of  IR.  Since  Whitney  uses  divided  differences  in  an  essential  way,  Ker¬ 
gin  looks  for  a  viable  generalization  of  the  divided  differences.  His  approach  retains  the 
univariate  choice  of  interpolating  from  7r^  at  an  arbitrary  (k  +  l)-set  T  in  IRd  and  deals 
with  the  many  more  degrees  of  freedom  available  from  nk  by  enforcing  certain  mean-value 
conditions.  These  conditions  are  that,  for  every  sufficiently  smooth  /,  the  interpolant  Pf 
should  be  in  7 r*  and,  for  every  r  <  k,  for  every  homogeneous  polynomial  q  of  degree  r,  and 
for  every  (r  +  l)-subset  W  of  T,  there  should  exist  a  point  in  the  convex  hull  of  W  at  which 
q(D)f  and  q(D)Pf  agree.  Here  (to  be  more  explicit)  q(D)  is  an  r-th  order  homogeneous 
constant  coefficient  differential  operator,  i.e.. 

q(D)  =  £  a(o)Z?° 

|U|  =  r 

for  certain  coefficients  a(o).  Surprisingly,  there  exists  exactly  one  linear  map  P  with 
these  properties.  This  linear  map  can  be  characterized  by  the  fact  that,  for  every  plane 
wave.  i.e..  every  function  /  of  the  form  /  :=  g  o  X  :  x  *  g[^x.  \)),  P  reduces  to  univariate 
interpolation  at  the  projected  point  set  (T,  A)  =  {(f.A;  :  t  £  T)  C  IR.  i.e.,  P(g  o  A)  - 
(F’xtq)  A,  with  P\t  univariate  interpolation  from  zr*. ( IR )  at  XT  =  T,  A).  In  other  words, 
for  a  plane  wave  /.  Pf(x)  is  the  value  at  x.  A)  €  IR  of  the  univariate  polynomial  of  degree 
at  most  k  which  matches  the  value  f{t)  at  (t.  A  .  all  /  e  T .  where  Hermite  interpolation  is 
used  in  case  of  coincident  points.  Full  understanding  of  this  process  (see  Micchelli.  1980) 
led  to  an  understanding  of  multivariate  B-splines.  of  which  more  anon. 
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It  seems  more  promising  to  give  up  on  polynomials  altogether  and  to  choose  the 
interpolating  function  space  5  to  depend  on  the  point  set  T  at  which  data  are  given.  The 
simplest  general  model  has  the  form 


-  t)c{t), 

t£T 

with  p  :  IRd  — *  IR  a  function  to  be  chosen  ‘suitably’.  Duchon’s  thin  plate  splines  (see, 
e.g.,  Meinguet,  1979)  use 

/  x  .  ,2 m-d  j  /n|x|,  n  even; 
r{z)  -  111  \  1,  n  odd, 

motivated  by  a  variational  argument,  while  Hardy's  multiquadrics  correspond  to  the 
choice 

p(x)  :=  yjl  +  jx|2. 

A  good  source  of  up-to-date  information  about  such  interpolation  methods  and,  in  par¬ 
ticular,  about  the  question  of  their  correctness,  is  the  recent  survey  article  of  Micchelli 
(1986). 


4.  The  B-form 

I  now  come  to  a  discussion  of  piecewise  polynomial  functions,  or  pp  functions  for 
short.  I  have  learned  from  the  people  in  Computer-Aided  Geometric  Design  that,  in  dealing 
with  smooth  pp  functions  on  some  triangulation,  it  is  usually  advantageous  to  write  the 
polynomial  pieces  in  barycentric-Bernstein-Bdzier  form,  or  B-form  for  short.  This 
form  relates  polynomials  to  a  given  simplex.  It  is  hard  to  appreciate  the  power  and  beauty 
of  this  form  because,  even  with  carefully  chosen  notation,  it  looks  forbidding  at  first  sight. 
Still,  I  want  to  point  out  its  structure  at  least. 

One  starts  with  a  (d+  l)-subset  V  of  IRd  in  general  position  and  considers  the  barycen- 
tric  coordinates  with  respect  to  it.  i.e..  the  Lagrange  polynomials  for  linear  interpolation 
at  V.  The  typical  Lagrange  polynomial  takes  the  value  1  at  the  vertex  r  and  vanishes 
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on  the  facet  spanned  by  V\v.  The  B-form  for  p  6  *k  employs  all  possible  products  of  k  of 
these  linear  polynomials.  Explicitly, 

p=:  B<>  c(a)  (4.1) 


£a(x)  :=  la|!|^(x)|a  =  |a|! 


I 


£(*) :=  (£»(*))  ..  (4-3) 

is  the  ( d  t  l)-vector  containing  the  barycentric  coordinates  of  x  with  respect  to  V. 

Note  that  the  vector  f  (x)  and  the  multi-index  a  appearing  here  are  conveniently  and 
appropriately  indexed  by  the  elements  of  V  (rather  than  by  the  numbers  1,2,...  ,d  +  1  or 
the  numbers  0, 1,. . .  ,d,  which  would  require  an  arbitrary  indexing  of  the  points  in  V). 

The  factor  |a|!  in  the  definition  of  the  Bernstein  basis  element  Ba  is  just  right  to 
make  (Ba)ia |_fc  a  partition  of  unity.  Indeed. 

Y.  Ba{x)  =  k\  y  n  iu*r(v)  =  k-i  Y  ^(x)ik  = 

|a|  =  fc  \a\  =  kveV  v£V 

using  the  Multinomial  Theorem  (see  (l.l))  and  the  fact  that  £v  =  1.  The  numerical 
analyst  will  delight  in  the  alternative  formulation  of  the  form, 

P(*)  =  <£(x),£)*c(0)  M) 

which  makes  use  of  the  shift  operator  E  given  by  the  rule 


m 

.*»'!• 


More  explicitly, 


ECic(a)  =  c(q  -  3). 


[z.E  ^  Vs(r)E,. 


.z.  E  c(a)  -  Y  2(i’)c(q  -  ct.) 


AV/ir.  rL*  \  A  -V.N  JVJV  kV  ^  kV  .  **VvVwV>V>*%. 
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with  tv  the  u-unitvector,  i.e,  ev(w)  =  Svw,w  G  V'.  This  form  provides  a  most  convenient 
starting  point  for  the  derivation  of  efficient  algorithms  for  the  evaluation  and  differentiation 
of  the  B-form.  For  details,  see,  e.g.,  Farin  (1985)  and  de  Boor  (1986). 

5.  Smooth  pp  functions 

The  B-form  is  well  suited  to  pp  work  since  its  typical  term  Ba  vanishes  a(v)-fold  on 
the  facet  spanned  by  V\v.  This  means  that  the  form  readily  provides  information  about 
the  behavior  of  p  at  all  the  bounding  faces  of  the  simplex  with  vertex  set  V’.  This  is  being 
increasingly  exploited  in  studying  the  algebraic  structure  of  the  space 

<A 

of  pp  functions  of  degree  <  k  on  a  given  triangulation  A  whose  pieces  join  together 
smoothly  to  provide  a  function  all  of  whose  derivatives  of  order  <  p  are  continuous. 

The  problems  being  studied  include;  the  dimension  of  such  a  space,  a  good  basis  for 
such  a  space,  and  the  approximation  power  of  such  a  space.  For  recent  results,  see  Chui 
and  his  co-workers,  and  Schumaker.  These  results  only  deal  with  d  =  2,  and,  even  for  this 
case,  we  know  relatively  little.  For  example,  despite  considerable  efforts,  we  still  do  not 
know  the  dimension  of  the  space  of  continuously  differentiable  piecewise  cubic  functions 
on  an  arbitrary  triangulation  in  the  plane.  While  we  do  know  that  this  dimension  depends 
on  the  quantitative  details  of  the  triangulation,  we  do  not  know  exactly  how. 

As  we  understand  these  problems  better  and  see  some  of  their  particular  difficulties, 
we  wonder  whether  7r£  A  is  really  the  right  space  to  study.  It  now  seems  that  it  might 
be  more  appropriate  to  seek  out  appropriate  subspaces,  e.g.,  the  subspace  spanned  by 
certain  compactly  supported  smooth  piecewise  polynomials  as  was  done  already  in  Finite 
Elements.  A  particularly  simple  model  is  provided  by  approximation  from  a  scale  of  pp 
functions. 


- 


:v:v* 


-VV1 
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6.  Approximation  from  a  scale 


Associate  with  a  given  function  space  5  the  scale  (5>,),  with 

Sh  :=  ohS,  ( of){x )  :=  f[x/h), 

and  define  the  approximation  order  of  5  to  be 

max  {r  :  V  smooth  /  dist (/,  5^)  -  0(/ir)}. 

This  order  may  well  be  0,  as  it  is  for  5  =  tt*.  But  if  S  contains  functions  whose  support 
has  diameter  6,  then  Sh  contains  functions  with  supports  of  diameter  /i6.  and,  for  such 
S,  one  might  hope  to  obtain  closer  approximations  from  Sh  as  h  — *  0.  Work  with  specific 
examples  has  suggested  the  following  conjectures  in  case  5  C  7rfc,A: 

Conjectures:  (i)  The  approximation  order  of  S  equals  the  approximation  order  of 

S|oc  :=  span  {tp  €  S  :  supp  <p  compact}. 

(ii)  S  has  approximation  order  >  1  iff  S  contains  a  local  partition  of  unity. 

(iii)  The  approximation  order  is  always  realized  by  a  good  quasi-interpolant. 

Here,  a  map  Q  into  5  is  a  good  quasi-interpolant  of  order  r  in  case  it  is  a  linear 
map  which  is  stable  in  the  sense  that,  for  any  /  and  any  jf  G, 

!(<?/)(*)!  <  const  sup{|/(y);  :  j;y  -  x||  <  R} 

with  const  and  R  <  oc  independent  of  /  or  x.  and  which  reproduces  polynomials  of 
degree  «-  r.  For  example,  if  4>  is  a  local  and  nonnegativc  partition  of  unity  in  5,  i.e., 
sup*,.  $  diam  suppip  <  oc,  <p  >  0  for  all  <p  i  4>,  and  -  1,  then 

/  ~  P  /(V) 

<t> 

is  a  good  quasi-interpolant  of  order  1  (provided,  e.g..  that  t ^  suppp  for  all  p  ►;  <l>). 
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This  abstract  model  can  be  completely  analysed  in  the  very  special  case  when  S  is 
spanned  by  the  integer  translates  of  one  function  < p,  i.e., 

5  :=  Sv  :=  spanf<p(*  —  j))  =  {  1P(*  -  3 Mi)  :  c(j)  6  IR}. 

For  this  case,  the  three  conjectures  are  verified;  in  particular,  Strang  and  Fix  (1973)  prove 
that  S  has  approximation  order  r  iff  it<r  C  S. 

Already  for  the  slightly  more  general  case  when  S  is  the  span  of  integer  translates  of 
several  compactly  supported  functions,  the  situation  becomes  more  complicated.  A  char¬ 
acterization  of  the  approximation  order  is  not  yet  known  for  this  case,  but  the  somewhat 
stronger  (and  practically  more  interesting)  concept  of  local  approximation  order  can  be 
characterized  very  simply  (de  Boor  and  Jia,  1985):  S  has  local  approximation  order 
r  iff  there  exists  x(j  e  Sjoc  such  that  has  approximation  order  r. 

For  the  general  case,  even  simple  questions  such  as  whether  a  pp  space  with  positive 
approximation  order  must  contain  a  compactly  supported  element  have  so  far  remained 
unanswered. 


7.  Multivariate  B-splines 

The  abstract  theory  of  approximation  from  a  scale  has  found  new  interest  recently 
because  of  the  advent  of  multivariate  B-splines.  These  were  introduced  in  1976  in  hopes 
that  they  would  perform  the  same  service  in  the  study  of  multivariate  smooth  pp  functions 
that  the  B-splines  of  Schoenberg  and  Curry  provided  so  nicely  for  the  theory  of  (univariate) 
splines. 

In  retrospect,  well,  in  any  case,  the  idea  is  simple  enough.  It  involves  a  body  B  C  IR" 

and  the  orthogonal  projector  P  :  IR"  — *  IRrf  :  x  — ►  (x(l) _ >x(d)).  The  map  P  is  used  to 

extend  a  function  ip  on  IR^  to  the  function 
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on  ail  of  IRn.  The  B-spline  Mb  is  defined  as  the  distribution  on  IRd  which  represents 
integration  over  B  of  the  extended  function.  In  formulae: 


Mb  '•  <p  •—>  /  <P  0  P 

Jb 


for  all  i p  6  Cq. 


Here,  C0  =  Co(IRd)  is  the  collection  of  all  continuous  functions  on  IRd  with  compact 
support.  If  PB  (the  projection  of  the  body)  is  d-dimensional,  this  can  also  be  written 


Mb{p)  =  /  MbP  =  /  dy4>{y )  / 
Jmd  Jpb  JL 


Br\P~  1  y 


showing  that  Msiy)  =  vo!n_<fi?  n  P~ly.  This  latter  formula  was  the  original  definition, 
motivated  by  a  geometric  characterization  of  the  (univariate)  B-spline  due  to  Curry  and 
Schoenberg  (1966)  and  illustrated,  for  n  -=  3,d  =  1,  in  Figure  7.1. 


Figure  7.1  The  quadratic  B-spline  as  the  “shadow”  of  a  3-simplex. 

The  value  of  the  B-spline  at  a  point  y  equals  the  (n  -  d)-dimensional  volume  of  the 
intersection  of  the  simplex  with  the  hyperplane  P~ly. 

It  is  immediate  that  Mg  has  compact  support.  Further,  if  B  is  polyhedral  with 
facets  {Bt},  and  if  z  6  fit”,  then  an  application  of  Stokes'  Formula  shows  that  the  direc¬ 
tional  derivative  of  Mb  along  Pz  is 

DpzMb  —  ~  ^  z.n,'M b,-  (8.1) 


mmmm 
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with  n,  the  outward  unit  normal  to  the  facet  Bt  and  Mb,  the  B-spline  that  is  the  “shadow” 
of  the  “body”  Bt.  Repeated  applications  of  this  differentiation  formula  show  that  all 
derivatives  of  Mb  of  order  n  -  d  +  1  must  vanish  identically  away  from  the  projections  of 
the  (d  -  l)-dimensional  faces  of  B.  Consequently, 

Mb  G 

with  A  the  partition  whose  partition  interfaces  are  the  projections  of  [d  -  l)-dimensional 
faces  of  B,  and  where  p  is  defined  by  the  condition  that  n-p-2  equal  the  largest  dimension 
of  a  face  of  B  projected  entirely  into  one  of  the  partition  interfaces.  Thus,  in  the  generic 
case,  we  have  p  =  n-d- 1,  which  is  as  large  as  it  can  possibly  be,  given  that  the  polynomial 
degree  of  Mb  is  n-d. 

This  surprising  smoothness  is  bought  at  a  price.  Since,  for  a  generic  partition  A, 
nk^A  does  not  contain  any  locally  supported  functions,  the  partition  for  Mb  must  be  quite 
special.  Figure  7.2  shows  such  a  partition  for  a  bivariate  quadratic  simplex  spline,  i.e.,  a 
B-spline  that  is  the  “shadow”  of  a  simplex. 


Figure  7.2  The  partition  for  a  bivariate  quadratic  simplex  spline 

Thus  we  cannot  expect  to  obtain  B-splines  for  every  partition.  At  best,  we  can  find  B- 
splines  whose  partition  refines  a  given  one.  For  the  case  of  simplex  splines,  such  a  collection 
of  B-splines  of  degree  k  can  be  constructed  rich  enough  to  provide  a  good  quasi-interpolant 
of  order  k  1.  There  are  even  stable  recurrence  relations,  found  by  Micchelli  (1980).  for 
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their  evaluation.  But  it  seems  that  their  use  is  computationally  quite  expensive  (see,  e.g., 
Grandine,  1986).  It  is  therefore  not  likely  that  simplex  splines  will  be  used  as  a  basis  for 
a  good  subspace  of  a  given  smooth  pp  space  of  functions.  Most  likely,  translates  of  a  fixed 
B-spline  will  find  practical  employment. 

There  is  a  bit  more  hope  for  the  box  splines,  i.e.,  the  multivariate  B-splines  associated 
with  the  n-cube  B  —  [0,1'".  For  their  definition  in  terms  of  (7.1),  one  would  allow  P  to 
be,  more  generally,  a  linear  map.  Then,  with  £,  :=  Ptx  the  image  of  the  t-th  unit  vector 
under  P,  the  box  spline  can  be  characterized  more  explicitly  by 


/  M(-|fi,  •••,£„)£>  =  (  <p(y2t,y{i))dy 


,  <p  €  Co- 


By  choosing  the  from  TLd  appropriately,  the  resulting  partition  can  be  made  to  conform 
to  a  regular  grid;  see  Figure  7.3  for  the  supports  of  two  very  well  known  box  splines,  the 
Courant  element  (d  =  2,n  =  3)  and  the  Zwart-Powell  element  (d  =  2 ,  n  =  4).  Further, 
their  evaluation  can  be  accomplished  by  subdivision. 


Figure  7.3  The  supports  of  a  linear  and  a  quadratic  bivariate  box  spline 

De  Boor  (1982)  gives  an  introduction  to  multivariate  B-splines.  Dahmen  and  Mic- 
chelli  (1984)  provide  a  survey  of  the  literature  available  by  1983  to  which  they  heavily 
contributed.  Hollig  (1986a)  gives  a  more  up-to-date  introduction,  and  Iioliig  (1986b)  sum¬ 
marizes  what  we  know  about  box  splines.  In  addition,  1  want  to  stress  the  beautiful,  but 


Multivariate  Approximation 


15 


more  theoretical,  developments  to  which  Dahmen  and  Micchelli  were  led  by  their  intensive 
study  of  box  splines  (see,  e.g.,  Dahmen  and  Micchelli,  1985,  1986). 

8.  Subdivision 

I  hate  to  finish  on  a  pessimistic  note.  I  therefore  bring  up  a  totally  different  approach 
to  the  generation  or  approximation  of  surfaces  which  comes  from  Computer-Aided  Design. 
I  think  that  this  technique  has  real  promise  for  the  generation  of  ‘smooth’  surfaces  which 
fit  to  given  points  in  3-space  of  more  or  less  arbitrary  combinatorial  structure.  It  is  at 
present  being  used  to  evaluate  linear  combinations  of  box  splines  (see,  e.g.,  Hollig,  1986a, 
for  the  relevant  references).  But  since  this  idea  has  not  yet  been  thoroughly  studied,  1  will 
discuss  it  only  in  its  original  context  of  curve  generation. 


Figure  8.1  The  steps  of  the  simple  subdivision  algorithm. 


Here  is  a  very  simple  version  of  subdivision,  which  generalizes  Chaikin’s  algorithm 
(Chaikin,  1974).  Start  off  with  points  a}  in  IRd.  where  j  runs  over  all  of  TL.  for  simplicity. 
Think  of  these  points  as  the  vertices  of  a  broken  line.  From  the  algorithm,  one  obtains  a 
refined  broken  line  in  two  steps.  In  the  first  step,  one  introduces  the  midpoints  between 
neighboring  vertices  as  new  vertices,  thus  roughly  doubling  the  number  of  vertices: 

62;  :=  a}.  l>2J  +  i  :=  (a}  a;* i)/2. 
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In  the  second  step,  one  obtains  each  vertex  of  the  refined  broken  line  as  an  average  of  three 
neighboring  vertices: 


c3  (3bj-\  +  abj  +  'y6;  +  i,  with  0  a  +  7  =  1. 


In  fact,  for  a  curve  of  higher  ‘smoothness’,  one  would  repeat  this  averaging  step  one  or 
more  times,  but  I  will  stick  with  this  simple  model.  Repetition  quickly  leads  to  a  broken 
line  which,  for  plotting  purposes,  is  indistinguishable  from  the  limiting  curve. 


Of  course,  it  is  not  at  all  clear  a  priori  that  there  is  a  limiting  curve,  though  that  is 
easily  proved  for  reasonable  choices  of  the  weights,  e..g,  for  a,  0,  7  >  0.  Nor  is  it  clear  just 
what  the  nature  of  that  limiting  curve  might  be.  Chaikin’s  algorithm  corresponds  to  the 
choice  0  =  0,  a  =  -7  =  1/2.  For  this  choice,  the  limiting  curve  is  a  parametric  quadratic 
spline  curve,  viz.  the  curve 

t  >-*  ^  Mz(t  -  j)aj 


with  Mz  a  quadratic  cardinal  B-spline  (i.e.,  a  B-spline  having  integer  knots).  For  the 
symmetric  choice  0  =  7  =  (1  -  a)/2,  the  limiting  curve  is  a  parametric  quadratic,  resp. 
cubic,  spline  curve  in  case  a  =  0,  resp.  1/2,  but  for  any  other  choice,  the  limiting  curve 
appears  to  be  something  unmentionable  in  standard  terms,  though  this  is  not  apparent 
from  the  curves  themselves;  see,  e.g.,  Figure  8.2. 


Figure  8.2  Curve  iterates  w  ith  a  =  1  A,  3  - 
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Since  the  limiting  curve  depends  linearly  on  the  initial  data  a;,  it  is  sufficient  for  the 
complete  analysis  of  the  process  to  consider  a  broken  line  with  just  one  actual  break,  e.g., 
the  data 

Uj  :=  j€7Z  (8.1) 

in  IR2.  The  various  curves  generated  will  differ  only  on  the  segment  between  the  two  points 
(—1,0)  and  (1, 1).  Some  of  the  iterates  as  well  as  the  limiting  curve  segment  are  shown  in 
Figure  8.3  for  the  particular  choice  a  =  1/4,  and  we  see  nothing  unusual. 


Figure  8.3  Iterates  and  limiting  curve  for  one-break  broken  line 


It  is  only  when  we  look  at  derivatives  that  we  realize  that  something  is  amiss.  For 
our  particular  curve,  the  first  component  is  just  t  t,  so  it  makes  sense  to  consider  the 
derivatives  of  the  second  component,  t  —*  y(f)  say,  as  an  indication  of  the  smoothness  of 
the  curve.  Since  we  do  not  have  a  formula  for  y(t),  we  cannot  compute  derivatives.  But 
since  the  iteration  is  so  simple,  we  can  compute  y(t)  for  as  many  binary  fractions  t  as  we 


18 


Carl  de  Boor 


care  to.  That  done,  we  can  then  compute  second  divided  differences,  and  these  tell  the 
story;  see  Figure  8.4:  The  second  derivative  of  y  appears  to  be  a  fractal. 

Traditional  Approximation  Theory  views  such  curves  with  horror.  There  is  even  the 
seemingly  practical  objection  that  it  would  be  difficult  to  machine  curves  with  such  a 
‘bad’  second  derivative.  But  I  do  not  know  of  experiments  which  have  established  such  a 
difficulty  nor  is  it  clear  a  priori  that  there  should  be  any  difficulty.  On  the  other  hand,  the 
generation  of  such  curves  to  plotting  or  machining  accuracy  is  so  swift,  and  their  flexibility 
so  great,  that  this  technique  is  worth  exploring  in  detail.  It  is  this  apparent  local  flexibility 
that  makes  subdivision  techniques  a  promising  tool  for  the  generation  of  shape-controlled 

surfaces  of  arbitrary  combinatorial  structure. 

*  * 


Figure  8.4  Wisconsin  Winter 

Second  divided  differences  it  —  h,t,t  -  h  y  plotted  against  f,  for 
h  -  2  B  (dots),  and  h  —  2'9  (stars) 

For  literature,  see  Catmull  and  Clark  (1978).  Doo  (1978).  and  Doo  And  Sabin  (1978).  I 
am  indebted  to  Professor  Wanner  for  the  surprising  references  to  work  by  de  Rham  (1947. 
1953,  1956,  1957,  1959)  who,  many  years  ago  and.  apparently,  for  pedagogical  reasons, 
investigated  a  related  subdivision  algorithm  for  curves.  Professor  Wanner  referred  to  it 
imaginatively  as  ‘the  woodcarver's  algorithm'. 
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